Mn Analysis

Check correlation between Mn responsive and microbe responsive genes in ime course

library(tidyverse)
library(readxl)
library(edgeR)

get the Mn data

Mn <- read_csv("../output/Mn_deficiency.DEGs.all.anno.csv") %>%
  select(AGI, logFC, FDR) %>%
  filter(FDR < 0.05)

Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  X1 = col_double(),
  AGI = col_character(),
  logFC = col_double(),
  logCPM = col_double(),
  LR = col_double(),
  PValue = col_double(),
  FDR = col_double(),
  symbol2 = col_character(),
  full_name = col_character()
)

Mn

Now get the microbe response data

load("../output/TimeCourseDGE_LRT.Rdata")

names(models)

[1] "day_time" "day"      "time"     "dge"      "lrt"      "DEG0.05" 

DEgenes <- models %>% 
  mutate(toptags = map(lrt, ~ topTags(.,n=Inf)), #pull the DE gene info out of each model
         toptags = map(toptags, ~ .$table), # get the table of DE info
         toptags = map(toptags, ~ rownames_to_column(., "transcript_ID"))) %>%
  unnest(toptags) %>%
  group_by(transcript_ID) %>%
  arrange(transcript_ID, FDR) %>% 
  select(day_time, transcript_ID, logFC, FDR)
DEgenes

Add Arabidopdis ID to each Brassica gene

First, get annotation

annotation <- read_csv("../../Brapa_microbes/Annotation/output/v3.0annotation/Brapa_V3.0_annotated.csv")

Missing column names filled in: 'X1' [1]Parsed with column specification:
cols(
  .default = col_double(),
  name = col_character(),
  chrom = col_character(),
  subject = col_character(),
  AGI = col_character(),
  At_symbol = col_character(),
  At_full_name = col_character(),
  At_gene_model_type = col_character(),
  At_short_description = col_character(),
  At_Curator_summary = col_character(),
  At_Computational_description = col_character()
)
See spec(...) for full column specifications.

annotation

now join the annotation

DEgenes <- annotation %>%
  select(name, AGI) %>%
  right_join(DEgenes, by=c("name" = "transcript_ID"))
DEgenes

now add in the microbe data

gene_compare <- DEgenes %>%
  inner_join(Mn, by = c("AGI"), suffix = c(".microbe", ".Mn"))
gene_compare

gene_compare <- gene_compare %>%
  group_by(day_time) %>%
  nest()
gene_compare

gene_compare <- gene_compare %>%
  mutate(cortest = map(data, ~ cor.test(.$logFC.microbe, .$logFC.Mn, method="spearman", exact=FALSE)))

gene_compare %>% mutate(glance = map(cortest, broom::glance)) %>% unnest(glance) %>%
  arrange(p.value) %>% select(day_time, estimate, p.value) 

Plot the most significant one

gene_compare %>%
  filter(day_time=="02-2_afternoon") %>%
  unnest(data) %>%
  ggplot(aes(x=logFC.Mn, y=logFC.microbe)) +
  geom_point(alpha=.2)